function cmask = RegularizeMRI(img,pixfrom,pixto,factor,filter_sd)
%% Regularize high-field MR image.   Calculate a 2D bias or the background noise filter across
% the main axis (1st dimension). The bias is taken from cross-sections of the
% volume without the sample and away from the ends.  The loss of SNR occurs across the third 
% dimension at positions furthest from the receiver coils. This backgound cross-section will be used
% to enhance values similar to the centre of the cross-section (moderate signal high SNR) 
% and reduce values with high signal and low SNR.   
%
% reg·u·lar·ize Verb: Make (something) regular. Establish (a hitherto
% temporary or provisional arrangement) on an official or correct basis.
%
% Note:  The tube containing the sample must be correctly aligned along the main axis. 
%
%  >> cmask = RegularizeMRI(img,pixfrom,pixto,factor)
%  >> corrected_img = MaskVolume(img,cmask);
% 
% Arguments
%  img       3D volume of raw image 
%  pixfrom   (optional) index of starting the slice average
%  pixto     (optional) index ending slice average
%             Without pixfrom/pixto the slice averaging method calculates points
%             using an estimate from a 2D histogram (across the second dimension).
%             The loss of SNR occurs across the third dimension.
%  factor    (optional) scaling factor of filter after normalisation
%  filter_sd (optional) standard deviation of the gaussian filter used on
%            the 2D mask
% Output
%  cmask     2D mask  in 3D volume that estimates the background bias field
%            along the main axis (1st dimension)
%
%   - Michael Eager,  (michael.eager@monash.edu)
%   - (c) 2012,  Monash Biomedical Imaging,  Monash University, Australia


%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

xmidpoint = floor(size(img,1)/2);
ymidpoint = floor(size(img,2)/2);

if nargin < 3
    disp('Estimating slices without kidney sample for RegularizeMRI mask.')
    disp('This method works for samples in a tube that is aligned along the main axis.')
    % get maxima along second dimension
    nmimg  =squeeze(max(img,[],2));
    % find maxima in midsection below 20000
    % 20000 is arbitrary and selected from example 1008.2.18.3.1.1.1
    % correction - now find the midpoint of the maxima 
    mean_along_xaxis_at_midsection = mean(nmimg(:,ymidpoint-50:ymidpoint+50),2);
    cross_point = ((max(mean_along_xaxis_at_midsection)-min(mean_along_xaxis_at_midsection)) /2 ) + min(mean_along_xaxis_at_midsection)
    indx = find(mean_along_xaxis_at_midsection  < cross_point);
    % only look at points in bottom half
    indx=indx(find(indx > xmidpoint));
    
    %Check size of indx
    if numel(indx) < 12
        pixfrom = indx(1);
        pixto = indx(end);
    else
        pixfrom = indx(4); % start five in from drop in maxima
        pixto  = indx(end-8); %ignore 8 slices from end
    end
    
    if pixfrom > pixto
        pixfrom = pixto;
    end
    
    disp(['RegularizeMRI slices without kidneys: ' num2str(pixfrom) ' to ' num2str(pixto) ])
   
    clear nmimg,indx
end
if nargin < 4
    factor=1;
end
if nargin < 5
    filter_sd=2.0;
end


nimg= NormaliseImage(img);
%Calculate background along main axis
bg =  squeeze(mean(nimg(pixfrom:pixto,:,:),1));
f = fspecial('gaussian', 10, filter_sd);
mask = imfilter(bg,f,'replicate');
mask = mask./max(mask(:));

% mean of centre region used as calibration point

cpoint = mean(mean(mask(xmidpoint-25:xmidpoint+25,ymidpoint-25:ymidpoint+25)))

% Final bias correction 2D mask
cmask = (abs(mask - cpoint ) .* factor);
